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ABSTRACT 


A new method of array modeling which will be used to predict 
the performance of low frequency active sonar arrays is being 
investigated. In support of this effort, a network representation of 
a spherical shell piezoelectric transducer was developed. The 
transducer was modeled using the finite-element code MARTSAM, 
from which a nodal description of the transducer was obtained. A 
procedure was developed to reduce and transform the nodal 
description of the transducer into a spherical harmonic description. 
The spherical harmonic description of the transducer was computed 
at two frequencies, 112.5 Hz and 1125.3 Hz, corresponding to values 
of ka of 0.1 and 1.0, respectively, where a is the radius of the 


sphere. 
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1. INTRODUCTION 


The objectives of the research detailed in this thesis are: 1) to 
produce a nodal description of a spherical shell, piezoelectric 
transducer using finite-element modeling, and 2) to develop a 
procedure to transform the nodal description of the transducer into 
a spherical harmonic description. 

This report is organized in three main sections. First, a new 
method to predict the performance of low frequency active arrays, 
which provided the motivation for this research, is described. Next, 
the finite-element modeling and description in terms of spherical 
harmonics of a particular spherical shell transducer are detailed. 
Numerical results are included for the transducer operating at 112.5 
Hz and 1125.3 Hz, corresponding to values of ka of 0.1 and 1.0, 
respectively, where a is the radius of the sphere. Lastly, the 
procedure to reduce and transform the nodal description of the 


transducer to a spherical harmonic description is described. 








li. MODELING OF LOW FREQUENCY ACTIVE ARRAYS 


A. BACKGROUND 

The trend toward operation of sonar surveillance systems at 
lower frequencies has necessitated larger, more dense arrays of 
transducers. This transition has prompted investigation of a new, 
potentially more efficient, method to predict array performance. 
This new method of array modeling combines a finite-element 
representation for each transducer with a mathematical 
representation of the acoustic field which is equivalent to the 
T-matrix method which has been applied to elastic scattering 
problems [Ref. 1]. This approach is explained in greater detail in - 
Appendix A. It has the feature that it accounts for the effects of 
‘multiple-scattering’, which can be quite significant in dense arrays 
(the distortion of the near-field radiation pattern of a transducer 
due to a nearby transducer, for example, is a manifestation of 
multiple-scattering). Applying this method to an array of 
transducers will ultimately yield a matrix which relates the 
outgoing acoustic waves to the input electrical excitation of each 
transducer. 

Figure 1 is a schematic diagram of a portion of an array. In 

this method, the modeling problem is divided into two regions. The 


term structure refers to the transducer plus some surrounding fluid 


out to an imaginary, spherical surface, S. The second region is the 
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fluid outside of S. The array may be composed of any number and/or 
variety of tranducers. Additionally, it may have any geometry. 

The goal is to develop an analytical representation of the 
structure from a finite-element model and to combine it with an 
analytical representation of the acoustic field in the fluid to predict 
the performance of an array of interacting transducers. By 
separately modeling the structure and the fluid, this approach will 
simplify the complex task of array modeling. Changing the geometry 
of an array will require recomputing only the fluid model. Likewise, 
maintaining the geometry and changing the type of transducer will 


require a new representation of only the sructure. 


B. FLUID REPRESENTATION 

The acoustic pressure field in the unbounded fluid is 
represented as an eigenfunction expansion in terms of spherical 
waves. It should be noted that all quantities are assumed to vary 
harmonically with time. This representation is convenient due to 
the spherical shape of the structures. Figure 2 is the radiation 
pattern which is produced by three pulsating spherical elastic shells 
of radius a, with one half wavelength separation and ka equal to 1.0. 
Each shell is driven with a uniform pressure amplitude on the inside 
of the surface of the strength indicated by the number inside. The 
curves shown represent the radiation pattern when the number of 
spherical harmonics retained in the expansion is varied. The curves 


which include the first 6, 9, and 11 harmonics are indistinguishable; 




















therefore, the acoustic field may be accurately described by 


retaining only about 6 spherical harmonics. ([Ref.2] 


C. STRUCTURE REPRESENTATION 

The structure must be spherical in shape; it may be a 
transducer with some surrounding fluid or simply a spherical 
transducer. The structure is first modeled using a finite-element 
code, which produces a nodal description of the structure. The nodal 
description is then reduced and transformed into a spherical 
harmonic description of the structure. This description can then be 
combined with the spherical harmonic description of the acoustic 


field on the surface of each sphere to predict the behavior of the 


array. 








lll. MODELING OF SPHERICAL SHELL TRANSDUCER 


A. FINITE ELEMENT MODEL 

A spherical shell piezoelectric transducer was modeled using 
the finite element code MARTSAM. Figure 3 is the finite element 
mesh used to model the transducer. The transducer is radially 
polarized Navy Type 1 ceramic with an outer radius of 7.620 
centimeters and thickness of 0.762 centimeters. It is symmetric 
about the y-axis. Based on this symmetry, it was only necessary to 
create a two-dimensional model of one half of the transducer. The 
structure was divided into 60 six-noded triangular elements. This 
partitioning scheme results in a mesh with 183 nodes. There are 
two mechanical degrees of freedom associated with each node 
except the six nodes along the y-axis. Boundary conditions were 
applied to set the motion of these nodes in the x-direction equal to 
zero. Consequently, this structure possesses a total of 360 
mechanical degrees of freedom. Additional boundary conditions 
included grounding the electrode on the outer surface of the sphere 
while maintaining a constant potential on the inner electrode. The 


structure, therefore, has one electrical degree of freedom. 


B. MARTSAM OUTPUT 
A modal analysis of the structure was performed by Dr. Michele 


McCollum of the Naval Research Laboratory, Orlando using the finite- 








Figure 3: Finite Element Mesh of Spherical Shell 
Transducer. 








element code MARTSAM. For this particular application, the code 
was used to generate (Kuu), (M), and Kyy which can be substituted 


into the following set of equations: 


a oS \@)- Bk (3.1) 


IV VV 


where 
« U and E are vectors that contain the nodal values of 
the displacement field and the applied forces, 
* V is the applied electrical potential, 
*« Q is the electrical charge on the structure, 
° ( Kuu) is the matrix which describes the effect that the 


stiffness at each node has on all the nodes, 

. (M) is the matrix, which describes the effect that the 
mass at each node has on all the nodes, 

* Kyy is the vector that contains the coupling coefficients 
that relate the mechanical and electrical degrees of 
freedom, 

*« Kyy is the capacitance of the transducer for zero 
displacement everywhere, 

¢@ is the angular frequency, 


¢ T means transpose, 














to completely describe the spherical shell transducer in terms of its 
nodes at a specified frequency. The matrix of the left side of 
equation (3.1) has dimension 361 by 361 for this structure. 

By means of matrix algebra, which is outlined in Chapter IV, 
this nodal description of the structure was reduced. The mechanical 
degrees of freedom associated with the nodes internal to the 
structure were eliminated as the external force on those nodes is 
zero. The displacements at the surface nodes were transformed 
from the Cartesian coordinate system to the polar coordinate 
system. The degrees of freedom associated with the parallel 
displacements were also eliminated as external fluid forces are 
applied normal to the surface. The dimensions of the reduced matrix 
are 62 by 62, 61 surface normal displacement degrees of freedom 


and 1 electrical degree of freedom. 


C. SPHERICAL HARMONIC DESCRIPTION OF STRUCTURE 
Following the procedure detailed in Chapter IV, the reduced 
nodal description of the transducer was transformed into a spherical 
harmonic description. Since the finite-element model is restricted 

to axisymmetric solutions for field quantities, the Legendre 
polynomials were used for the spherical harmonic functions. Based 
on the results obtained by Canright and Scandrett [Ref. 2], the first 7 
Legendre polynomials were used. The spherical harmonic description 
of the structure is an 8 by 8 matrix. Figures 4 and 5 are the 


electrodynamic matrices obtained for the frequencies 112.5 and 








1125.3 Hz, respectively, corresponding to values of ka of 0.1 and 1.0, 


respectively, where a is the radius of the structure. These 
matrices, (Ds), are used in the following set of equations: 


(psn) (YY) = es! (3.2) 


where USPh and oSPh are vectors that contain the values of the 
displacement field and stress field in terms of Legendre polynomials 
(n=0,1,2...6), to describe the transducer. Most of the elements in 
these matrices are on the order 108-1019. These numbers 
represent mechanical interactions while the elements in the last 
row and column (on the order of 10-5 - 10) represent the mechanical 
and electrical coupling. The element in the last row and column is 


the (surface only) blocked capacitance. 
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IV. TRANSFORMATION PROCEDURE 


A. REDUCED NODAL DESCRIPTION 


Matrix equation (3.1), for simplicity, can be represented by 
(Dnod) Ynod = Frod, (4.1) 


where 
° ( D>) is the square dynamical matrix, 
¢ Unod js the vector containing the displacement field 
and the electrical potential, and 
¢ Fnod js the vector containing the applied force field 


and the electrical charge. 


The dimensions of the dynamical matrix are quite large for even the 
simplest transducer. In modeling the fluid-loaded structure, 
however, the size of this description can be greatly reduced by 
applying matrix algebra as the following procedure describes. In 
performing this reduction, the accuracy of the description is not 
diminished since the fluid forces act only at the surface bounding 
the structure. 
1. Forces and Internal Nodes 
Fnod js initially a column vector whose entries are the force 


applied to each node listed by ascending node number and direction 


14 





followed by the electrical charge. However, the applied force is 
zero at all the nodes that are internal to the structure as those 
nodes are inaccessible to the surrounding fluid. [4 is rearranged so 
that all the zero entries come first. When F%¢ , the vector on the 
right hand side of the equation, is reordered, the rows of the 
dynamical matrix on the left side of the equation must be rearranged 
accordingly. U"°d must also be reordered in the same manner as Fnod 
which necessitates a corresponding rearranging of the columns of 
the dynamical matrix. 
2. Transformation of Coordinate System 
MARTSAM describes the transducer in terms of Cartesian 


coordinates. A transformation matrix, (T), can be defined for each 


surface node such that it performs the following operations: 


(us | (T » (er) and 


(Ex) = (ny) (er), Pe 


where, respectively, 
¢ U, and Uy are the displacements in the x- and y-directions, 


¢ Ex and Fy are the applied forces in the x- and y-directions, 
y 
* Unerp and Upar are the displacements in the directions 


perpendicular and parallel to the surface, and 
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¢ Eperp and Fpar are the applied forces in the directions 


perpendicular and parallel to the surface. 


For the mesh shown in Figure 3, for a single element, 
sin@  ~—cosé xX -Y 
(T)= ( cose — sino )- we ( Y XxX ), 48) 


where 
¢ @ is the angle each surface node makes with the positive 
y-axis, 
¢« X and Y are the Cartesian coordinates of each node, and 
¢ R is the radially distance of each node from the center of 
the sphere. 


Substituting equations (4.2) into equation (4.1) yields 


(Dred) (T) Upolar = (T) Epolar, (4.4) 


where the vectors U and F are now expressed in terms of polar 


coordinates. Multiplying each side of equation (4.4) by the inverse of 
(T) gives 


(1)-+ (De) (7) Upot -Epoa, 8) 
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where (T)-1 ( Dd) (T) represents the dynamical matrix in terms of 


nodes and the polar coordinate system. 
3. Forces and Parallel Direction 
The elements of FP°lar are arranged as follows: the forces 
applied to the internal nodes, the forces applied in the directions 
perpendicular to the transducer at the surface nodes and parallel to 
the transducer at the surface nodes, and the electrical charge. 
However, since only the forces applied normal to the transducer are 
non-zero, Fpar is the null vector. Using the procedure described in 
section 1, FPolar and Ypolar must be reordered such that Fpar and Upar 
precede Foerp and Uperp. The rows and columns of the dynamical 
matrix must be rearranged accordingly. 
4. Reduction of Dynamical Matrix 
The procedures in sections 1, 2 and 3 have arranged the 
dynamical matrix into the proper form to be reduced using matrix 


algebra. The matrix equation is now in the form 


(UL) Lay ) (upc) (Eber): (46) 


where 
. ( UL),(UR), (LL), (LR) are subdivisions of the 
dynamical! matrix, 
* Uo is the displacements at the internal nodes and in the 


direction parallel to the surface nodes, 


17 








* Unerp is the normal displacements and the electrical 


potential, 
* Q is the null vector, and 
* Foerp are the forces applied normal to the surface and the 


electrical charge. 

The corresponding simultaneous equations are 

(UL) Uo + (UR) Uperp = 0 , and (4.7) 

(LL) Uo + (LR) Uperp = Eperp - (4.8) 
Solving equation (4.7) for Uo yields 

Uo =- (UL)-1 (UR) Uperp - (4.9) 
Substituting equation (4.9) into equation (4.8) gives 

( (UR) - (LL) (UL): 1(UR) ) Userp = Eperp, (4.10) 


where ( (LR) - (LL) (UL)- 1(uR) ) is the reduced dynamical matrix 


describing nodal interactions. This matrix is square, and its 
dimension is the number of surface nodes plus the electrical degree 


of freedom. The computer code, which computes the reduced 
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dynamical matrix for the spherical shell transducer described in 
Chapter Ill, is listed in Appendix B [Ref. 3:p. 38]. 


B. TRANSFORMATION TO SPHERICAL HARMONIC DESCRIPTION 
The reduced nodal description of the transducer given by 


equation (4.10) can be written as 


(4.11) 


[feet ES JO) -s') 


KnodT Ce Q 


where 

- ( D8") is the reduced (Ku) - @2 (M) matrix, 

* Knod is the reduced vector containing the coupling 
coefficients that relate the mechanical and electrical 
degrees of freedom, 

* Ces is the (surface only) blocked capacitance, 

¢ UNod jis the vector containing normal displacements, 

« V is the electrica! potentiai, 


* Fnod js the vector containing the forces applied normal to 
the surface, 


* Q is the electrical charge, and 


¢ T means transpose. 


A spherical harmonic description of the transducer can be obtained, 
in the form 
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(Dseh) ysph = gsph | (4.12) 


by performing the following tranformations. 
1. Transformation of Yd to Usph 
The surface normal displacement distribution can be written 


as 


N 
U(e) = LUA NA(@) , (4.13) 
n 


where UR = U(® = On) , and N(6) is the finite-element interpolation 
function where Np(6) = 5(6 = 8n) [Ref. 4]. Since the finite-element 
model created to represent the trarscucer in this investigation is 
restricted to axisymmetric soiutions for all field quantities, the 
Legendre polynomials were chosen for the set of harmonic functions. 
The surface normal displacement distribution can then be expressed 


approximately by a truncated series as 


M 
U(8) = Zui Pm(cosé) , (4.14) 


T 
where Unp” = (1/Cm) {U(@) Pm(cose) sine de and Cm = |(Pm)2 sine dS, 
0 


the usual normalization constant associated with Pm. For the 
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Legendre polynomials, Cm = (2Il+1)/2 where | is the order. From 
equations (4.13) and (4.14), it follows that 


M 
h 
une? = DLUSP" B. (cose=cosén) . (4.15) 
m 


Equation (4.15) expressed in matrix form is 


Yrod = (P) sph, (4.16) 


where the matrix elements Pom = Pm(cos@=cosé6,) . 
2. Transformation of F"°d to gsph 
The self-consistent force at node n due to a distributed 


stress o is given by 


Frned = [oNndS, (4.17) 


where the integration is performed over the surface of the 
Structure. The distributed stress can be expanded similarly to the 


displacement as 


M 
0(8) = OFF" P.,(cose) . (4.18) 
m 
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The dynamical equation describing the transducer in terms of g§Ph 
follows from the equation for the coefficients of the expansion of 


the stress field in the Legendre polynomials: 


on” = (1/Cm) Jo Pm sine dS. (4.19) 
The expansion, 
N 
Pm(®) = >Pnm Nn, (4.20) 
n 


where Pram = Pm(8n), can be made. Substitution of equation (4.20) 


into equation (4.19) gives 


N 
on” = (1/Cm) D,  (foNndS) Pan. (4.21) 
n 


Using equation (4.17), equation (4.21) becomes 
N 
oP” = (4/Cm) UFR! Pron. (4.22) 
n 


In matrix form, equation (4.22) can be written 


gsph = (C)-1 (P)T Fnod , (4.23) 
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where 
+ (C)= diag (Cn). 
° (P) nm = Pm(cos@n) , and 
> FA = foNndS. 


Equation (4.23) is the desired expression relating g§P4 to Fnod, 


Matrix equation (4.11) can be expressed as 
(D884) ynod + Knod y = Fned | and (4.24) 
Knod T ynod + Cen V=Q. (4.25) 
Substitution of equation (4.24) into equation (4.23) gives 
gsph = (C)-1 (P)T (DEE%) ynody (C)-1 (P)T Kno. (4.26) 


Therefore, by substituting equation (4.16), equation (4.26) reduces 


to 
(DER") ysph 4 KsPh Vv = sph, (4.27) 


where (Dtb") = (c)-1 (P)T (DU8°) (P) and Keen = (C)-1 (P)T Krod. 


It follows from equations (4.16) and (4.25) that 


Krod T (P) UsPh + Cen V=Q. (4.28) 








Hence the reduced dynamical matrix equation in terms of 


Legendre polynomials coefficient degrees of freedom is 


ciety (OF (EH) (90) 


The computer code, which computes the transformation from the 
reduced nodal dynamical matrix to the spherical harmonic matrix, is 
listed in Appendix C (Ref. 3:p. 38]. 
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V. CONCLUSIONS AND RECOMMENDATIONS 


This thesis represents the first step in an investigation which, 
when successfully completed, will significantly improve the Navy's 
ability to predict the performance of dense sonar arrays. A network 
(matrix) representation of a spherical shell transducer has been 
calculated in terms of spherical harmonics. The procedure for 
reducing and transforming a nodal description of a transducer into a 
spherical harmonic description has been documented. 

For this thesis, a two-dimensional mode! of a spherical shell 
transducer was created and a procedure was developed to reduced 
and transform its nodal description into a spherical harmonic one, it 
is recommended that a three-dimensional model be created and a 
similar transformation procedure be developed. Since it is desirable 
that w remain a free parameter when describing a transducer, it is 
suggested that the procedure detailed in Chapter IV be rederived 
following Benthien's method to obtain the dynamical matrix as a 
function of » [Ref. 5]. It is also recommended that experimental data 


be obtained to verify the matrix representation developed in this 


thesis. 











APPENDIX A: 
EXCERPTS FROM PROPOSAL FOR LFA ARRAY MODELING 


APPLICATION OF THE T-MATRIX METHOD TO LOW-FREQUENCY 
ACTIVE ARRAY PERFORMANCE PREDICTION 


TECHNICAL PROPOSAL 


A. Introduction and Related Work 


The direction of active sonar surveillance systems is toward lower frequencies, requiring 
arrays of large, high power transducers. The successful design and operation of such arrays 
requires the ability to reliably predict their performance. 


Array performance prediction is a coupled structure-fluid problem [1]. For the analysis of 
complex structures, such as a sonar transducer, the finite-element method (FEM) is preeminent 
(2]. Many FEM computer codes exist, some of which include piezoelectric elements for the 
analysis of piezoelectric transducers, e.g. MARTSAM, ATILA. For application to coupled 
structure-fluid problems, the major drawback of the finite-element method is the difficulty of 
modeling the infinite fluid. So-called boundary elements, derived from a Helmholtz Integral 
formulation of the exterior radiation problem, have been introduced into the FEM to terminate 
the fluid mesh [3,4,5]. These elements typically impose an outgoing radiation impedance 
condition on the adjacent fluid element. The frequency dependence of the radiation impedance 
is commonly approximated by an interpolation between the low- and high-frequency 
asymptotic limits [5]. The application of a combined FEM-boundary integral method to array : 
performance prediction has not appeared in the literature. A full Helmholtz Integral formulation 
can be applied at the bounding surface [2,3], but this technique is limited to small problems by 
the number of degrees of freedom, and so is not practical to model an entire array. Hence there 
is a need for economical methods for array performance prediction. This proposal addresses 
that need. 


The method we propose is formulated in terms of coupled networks, one for each radiating 
element and one for the acoustic field. It is equivalent to the T-matrix method which has been 
applied to scattering problems [7,8,9,10], in that ultimately a matrix can be derived which 
relates the outgoing acoustic wave to the input electrical excitation of each transducer (the 
equivalent reciprocal problem is to relate the.electrical output of each transducer to an incident 
wave). One could in fact consider the proposed research the application of the T-matrix 
method to sonar array performance prediction. As in the T-matrix method, we adopt the idea 
of representing the acoustic field as an eigenfunction expansion and solving for the coefficients 
self-consistently. No restriction is placed on the arrangement of the radiators, and so multiple- 
scattering of all orders is rigorously included. In addition, we plan to explicitly take into 
account the dynamical properties of the radiators. These are to be found for a rea’ transducer 
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using a finite-element computer code [2,13]. 


B. Proposed Effort 
Goal 


The goal of the proposed research is to produce the most economical yet complete 
description possible of sonar array performance, with specific application to low frequency 
active arrays. The idea is to combine an analytical representation of a single element (derived 
from a finite-element model for a real transducer) with an analytical representation of the 
acoustic field to predict the performance of an array of interacting transducers. Especially for 
the usual case of an array of identical elements, since the dynamic behavior of only a single 
element need be computed, and since this computation need be done only once, regardless of 
array geometry, a variety of operating frequencies and array geometries may be analyzed most 
economically. 


Method 


A schematic diagram of a portion of an array of sonar transducers is shown below. 


ransducey 
Gronsducs} 


"Structure"-"fluid" boundary 





Transducey } 


"Structure' 













"Fluid" 


Figure 1. Schematic of a portion of an array of sonar transducers. 
The transducers need not be identical, nor in any particular arrangement. There may or may 
not be an acoustic wave incident from a source external to the array. The closed surface labeled 


S indicates the boundary between "structure" and "fluid". S may coincide with the physical 
ooundary of a transducer, or it may include some surrounding fluid. 
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There are several reasons to choose S to lie within the fluid surrounding a transducer. The i 
most important is that, if the "structure"-"fluid" boundary coincides with a constant-coordinate 
surface of a coordinate system in which the wave equation is separable, then the eigenfunctions 
of the wave equation in that coordinate system are a particularly convenient choice as the basis 
for expanding both the surface normal velocity on S and the acoustic field within the "fluid". 
This makes it much easier, for example, to compute the self- and mutual radiation impedances. 
Also, by including some surrounding fluid as part of the "structure", the burden of carrying 
many higher-order (higher spatial frequency) terms in an eigenfunction expansion of the 
acoustic field, which might be required to describe the local flow field near a real transducer, 
would be relieved. These degrees of freedom would be considered as internal to the structure, 
with the result that far fewer degrees of freedom would be required to describe the acoustic 
field. This is important because only the description of the acoustic field within the "fluid" 
needs to be recalculated for a change in array geometry. 


Conversely, by including some surrounding fluid, the description of the acoustic field 
would be unchanged by a change in the physical structure, for example, by exchanging one 
type of transducer for another (e.g. flextensional for hydroacoustic). The "structure" could 
even be composed of more than one physical transducer. This would be useful, for example, 
to analyze the performance of an array composed of groups of closely-spaced transducers. 


Including some surrounding fluid as part of the "structure" raises the possibility of 
introducing spurious structural resonances, which may cause problems for a numerical 
computation [10]. The consequences of this will have to be investigated. 


It is convenient to discuss the array performance problem further in terms of a coupled 
network representation. The diagram below depicts the coupling between "structure" and 
"fluid" for one surface normal velocity degree-of-freedom. For simplicity, the transducer is 
considered to be reciprocal. 
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Surface normal velocity 
degree-of-freedom n 
Transducer i Acoustic field 


in 





Figure 2. Coupled reciprocal network representation of structure-fluid 
interaction for one surface normal velocity degree-of-freedom. 


The network equations representing transducer i and the acoustic field are 


e= Zgiij - DTpuix, (transducer i) 
k 
fn= Twi - 2 Zr inktike 
fin=  DiApet % Ze injm4jm- (acoustic field) 
j.m 


The subscripts i and j refer to particular transducers, and the subscripts k, n and m refer to 
particular surface normal velocity degrees-of-freedom (DOF) on a “structure"-“fluid” 
boundary, assumed to be a constant-coordinate surface of a coordinate system in which the 
wave equation is separable. The meaning of the remaining symbols is 


e; = the voltage across the electrical terminals of transducer i, 
i; = the current through the electrical terminals of transducer i, 
fin = the modal force amplitude of transducer i, DOF n, 
Un = the modal surface normal velocity amplitude of transducer i, DOF n, 
Zz, = the blocked electrical impedance of transducer i, 
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Zm,inn = the open-circuit self mechanical impedance of transducer i, DOF n, 
Zm,inm = the open-circuit mutual mechanical impedance between transducer i, DOF n, and 
DOF m 
T,,, = the transduction coefficient of transducer i, DOF n, 
Z,.inin = the self radiation impedance of transducer i, DOF n, 
Z;,injm = the mutual radiation impedance between transducer i, DOF n, and transducer j, 
DOF m, 
Py= the incident free-field pressure, if any (usually assumed to be a plane wave), 
A; = the surface area of the "structure"-"fluid" boundary of transducer i, 
Din = the diffraction constant of transducer i, DOF n (the diffraction constant equals the 
ratio of the blocked modal force amplitude per unit area to the incident free-field 
pressure). 


The acoustic field is to be represented as the superposition of the acoustic field due to each 
radiator and an optional incident field. The surface normal velocity and the acoustic field of 
each radiating element are to be represented as expansions in the free-space eigenfunctions of 
the wave equation. *. ddition theorem [8,11,12] is to be used to express the field due to one 
radiator at the surface of another. The Z, injm and Dip are to be found in terms of the 
eigenfunction expansions by applying the boundary condition at the surface of each radiator 
that Ujm=5jOmn aNd Ujn=0, respectively. No restriction is placed on the arrangement of 
the radiators, and so multiple-scattering of all orders is rigorously included. A self-consistent 
solution for the coefficients of the eigenfunction expansions for each radiator is to be found by 
applying an impedance-matching boundary condition at the surface of each radiator. The 
mechanical impedance of each radiator is to be represented in terms of the modal parameters of 
its in-vacuo eigenfunctions. These are to be found for a real transducer using a finite-element 
computer code [2,13]. 
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APPENDIX B: 


PROGRAM REDDYNMTRX: computes the reduced nodal dynamical matrix 
for LT McLean's radially polarized piezoelectric spherical shell 


real omega,omega2,pi,freq,rce,pint( 361,361) 

real dyn(361,361),kuu(360, 360) ,mas( 360,360) 7 
teal r(361,361),9(361,361) 

real a(299,299),b(299,62),c(62,299) ,d(62,62) 

real nn(62,62),sum 

real p(62,299),pp(62,62),prod( 361,361) 

ceal s(361,361),mat( 361,361) ,y(299,299),indx(299) 
real t(361,361),tt( 361,361) 

real phtoh(8,86) 

real ccos(59),ssin(59),c2(59) 


OLD DATA FILES 
open(unit=5,file=’transform.sph’ ,status=’old’ ) 
open(unite-8, file="kuu.sph',status#'old’ ) 
open(unit=9,file='mas.sph’ ,status='old’ ) 
open(unit=10,file=’kue.sph’,status=’old’) 
open(unit=11,file=’kuel.sph’,status='old’ } 


NEW FILE FOR STORING THE REDUCED DYNAMICAL MATRIX 
open(unit=4,file='reddynmtrx.dat’,status=’new’ , form=’unformatted’ ) 


pi=acos(-1.0) 
ASSEMBLE THE FULL DYNAMICAL MATRIX 


Read in Kuu(360, 360) 
do i=1,360 
do j=1,360,5 
read(8,*)(kuu(i,k), k=j,j+4) 
end do 
end do . 


Read in Muu(360, 360) 
do i=1,360 
do j=1,360,5 : 
read(9,*)(mas(i,k), k=j,j+4) . 
end do 
end do 


Read Kue(360,1) into the 361th column of dyn( 361,361) 
do i=1,360 

read(10,*)dyn(i,361) 
end do 


Read Kue transpose (1,360) into the 361th row of dyn(361,361) 
do j=1,360 

cead(11,*)dyn(361,j) 
end do 


dyn(361,361)=capacitance for zero displacement everywhere 
dyn( 361, 361)=1.8515495E-08 


ASSEMBLE K-w**2M MATRIX 
write(*,*)’Enter frequency in Hz:' 
read(*,*)freq 
omega=(2.0*pi*freq) 
ome ga2=omega**2 
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do i=1,360 
do j=1,360 
dyn(i,j})=kuu(i,j)-omega2*mas(i,j) 
end do 
end do 


TEST SYMMETRY OF ASSEMBLED FULL DYNAMICAL MATRIX 
nerc=0 
do i=1,361 
do j=#i,361 
if (dyn(j,i).ne.dyn(i,j)) then 
write (*,*)'Symmetry error in full dyn. matrix at (i,j)=",i,j 
nerr=nerr+t 
end if 
end do 
end do 
write (*,*) ‘Number of symmetry errors in full dynamical matrix =’,nerr 
if (nerr.gt.0) then 
STOP 
end if 


REARRANGE ROWS SO ZERO (INTERNAL) FORCES COME FIRST 

do j=1,361 

r(1, j)=dyn(1, j) 

£(2,j})=dyn(2,3) 

(241, j)=dyn(3,j) 

r(3,j)=dyn(4,3) 

t(4,4)=dyn(5,j) 

£r(5,j)=dyn(6,}3) 

r(6,3)=dyn(7,3j) 

r(242,j})=dyn(8,j) 

c(243,j)=dyn(9,3) 

r(7,j)=dyn(10,j) 

r(8,})=dyn(11,j) 

r(9,j)=dyn(12,3) 

£(20, jdedyn(13, ) 

£(244,4)=dyn(14 j 

c(245,4)=dyn(15,j 

r(dl, j)edyn(16, ) 

r(12,j)edyn(17, j) 

£(13,j)=dyn(18,3j) 

r(14,j)=dyn(19,35) 

r(246,j)=dyn( 20,3) 

£(247,j)=dyn(21,3) 

r(15,4)=dyn(22,j) 

6(16, j)=dyn(23,j) 

r(17,j)=dyn(24,3) 

r(18,4)=dyn(25,j) 

£(248, j)}=dyn(26,j) 

r(249,4)=dyn(27,j) 

r(19,j)=dyn(28,j) 

£(20,j})=dyn(29,j) 

£(21,})=dyn( 30, 3) 

c(22,j)edyn( 31,4) 

(250, j)=dyn(32,j 

r(251,})=dyn(33,3 

r(23,j)=dyn( 34,3) 
) 
) 
) 


) 
) 


J 
Se ae die ae 
(25, j)=dyn(36,j 
(26,4) =dyn( 37,3 
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r(252,j)=dyn(38,j) 
v (253, })=dyn(39,j) 
c(27, j)=dyn(40,j) 
c(28, j)=dyn(41,j) 
r(29,j)=dyn(42, 3) 
£(30,4)=dyn(43,4) 
£(254,j)=dyn(44,j) 
c(255,})=dyn( 45,3) 
r(31,j)=dyn( 46, 4) 
r(32,j)=dyn( 47,3) 
£(33,j)=dyn(48,j) 
t(34,j)edyn( 49, 3) 
c(256,j)=dyn(50,j) 
r(257, })=dyn(51,4) 
r(35,43)=dyn(52,4) 
r(36,j)=dyn(53,j) 
1(37,j)=dyn(54, 3) 
c(38,j)=dyn(55, 4) 
r(258,j)=dyn(56,3) 
r(259, j)=dyn(57, 3) 
r(39,j)=dyn(58,j) 
r(40,j)=dyn(59,3) 
c(41,j)=dyn(60,j3) 
r(42,3)=dyn(61, 3) 
(260, 4)=dyn(62,}) 
r(261,4)=dyn(63,3) 
r(43, 3) =dyn(64,4) 
r(44,j)=dyn(65,j) 
r(45,4)=dyn(66,j) 
£(46,4)=dyn(67,j) 
r(262,4)=dyn(68,j) 
r(263,4)=dyn(69,}) 
r(47,j)=dyn(70,3) 
c(48,j})=dyn(71,3) 
r(49,j)=dyn(72,j) 
r(50,j)=dyn(73,3) 
r(264,4)=dyn(74,j) 
r(265,j))=dyn(75,j) 
r(51,j)=dyn(76,j) 
c(52,j)=#dyn(77,4) 
r(53,j)=dyn(78,j) 
r(54,4)=dyn(79,3) 
c(266,j)=dyn(80, 3) 
r(267,4)=dyn(81,j) 
r(55,j)=dyn(82,3) 
r(56,3)=dyn(83,j) 
r(57,j)=dyn(84,j4) 
r(58,5)=dyn(85,4) 
r(268,4)=dyn( 86,3) 
(269, j}=dyn(87,j) 
r(59,j)=dyn(88,j4) 
£(60,j)=dyn(89,4) 
c(61,j})=dyn(90,j) 
c(62;j)=dyn(91,3) 
r(270,4)=dyn(92, 4) 
r(271,4)=dyn(93,3) 
6(63,j)=dyn(94,j) 
r(64,4)=dyn(95,j) 
(65,j)=dyn(96,j) 
£(66,j)=dyn(97,j3) 
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£(272,4)=dyn(98,j 
(273, j)=dyn(99,j 
r(67,4)=dyn(100, 5 
c(68,4)edyn(1i01,j 
r(69,j)=dyn(102,j 
r(70,3)=dyn(1i03, j 
r(274,j)=dyn(104, 
£(275,4)=dyn(105, 
r(71,j)=dyn(106, j 
(72, j)=dyn(107,j 


) 
) 
) 
) 
) 
) 
i) 
j 
) 
}) 
c(73,j)=dyn(108, 3) 
}) 
J 
j 

) 
) 
) 
) 
j 


j) 


r(74,4j)=dyn(109,j3 
c(276,4)=dyn(110,4) 
(277, j)=dyn(111,j) 


£(77,j)edyn(114, 


r(278,j)=dyn(116 
£(279,4)=dyn(117,j 
r(79,j)=dyn(118,j) 
r(80,j})=dyn(119,j) 
(81, j)=dyn(120,3) 
r(82,})=dyn(121, 4) 
(280, 4)=dyn(122, 3) 
r(281, })=dyn(123,j) 
r(83,4)=dyn(124,j) 
r(64,j)=dyn(125,j) 
£(65,j)=dyn(126,j) 
r(86, 3)=dyn(127,j) 
c(282, 53) =dyn(128,j) 
r(283,3)=dyn(129,j) 
r(87,4)=dyn(130,j) 
r(88,j)=dyn(131,3) 
(89, 3) =dyn(132,j) 
r(90,j)=dyn(133,3) 
(284, 4)=dyn(134,j) 
r(285,j)=dyn(135,j) 
r(91, 4) =dyn(136,3) 
r(92,})=dyn(137,j) 
r(93,j)=dyn(138, 4) 
r(94,j})=dyn(139,j) 
r(286,j)=dyn(140,j) 
r(287,j)=dyn(141,j) 
£(95,j)edyn(142,4) 
£(96,j)=dyn(143,j4) 
£(97,j)edyn(144,j5) 
c(98,j)edyn(145, 3) 
r(288,j)=dyn(146,j) 
r(289, j)=dyn(147,j) 
£(99,j)=dyn(148,j) 
r(100, j)=dyn(149,j) 
r(101,4)=dyn(150,j3) 
6(102, j)=dyn(151,j) 
r(290,4}=dyn(152,j) 
(291, 4)=dyn(153,j) 
r(103,j)=dyn(154,j) 
r(104,4)=dyn(155,j) 
£(105, 4) =dyn(156,j) 
£(106,4)=dyn(157,j) 


) 
) 
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r(292,4)=dyn(158,3) 
r(293,j)=dyn(159,}) 
r(107,4)=dyn(160, 3) 
r(108, j)=dyn(161,j) 
5(109, })=dyn(162, 3) 
c(110,3)=dyn(163,j) 
r(294,4)=dyn(164,j) 
r(295,4)=dyn(165,j) 
r(111,3)=dyn(166,j) 
r(112,})=dyn(167,j) 
r(113,4)}=dyn(168, 3) 
r(114,j)=dyn(169,3) 
r(296,4)=-dyn(170,j) 
£(297,4)=dyn(171,4) 
r(115,4)=dyn(172, 4) 
r(116,4)=dyn(173,j) 
(117, j)=dyn(174,j) 
£(118, j)=dyn(175, 4) 
c(298,4)=dyn(176, 3) 
r(299, j)=dyn(177,4) 
(119, 4)=dyn(178, 3) 
£(120,4)=dyn(179, 3) 
r(121,4)=dyn(180, 3) 
r(122,})=dyn(181, 4) 
r(300,j)=dyn(162, 3) 
r(301, })=dyn(183,j) 
r(123,4)=dyn(184,j) 
r(124,4)=dyn(185, 4) 
r(125, j)=-dyn(186, 3) 
£(126,4)=dyn(187,j) 
r(302,4)=dyn(188, 4) 
r(303,4)=dyn(189, 3) 
(127, 4)=dyn(190,j) 
c(126,4)=dyn(191,4) 
c(129, 5) =dyn(192,3) 
r(130,4)=dyn(193,j) 
r(304,j)=dyn(194,j) 
(305, j)=dyn(195, 3) 
r(131,4))=dyn(196,3) 
(132, 4)=dyn(197,j) 
r(133,j)=dyn(198,j) 
r(134,4)=dyn(199,j) 
r(306,4)=dyn(200,j) 
c(307,4)=dyn( 201, 3) 
6(135,4)=dyn( 202, 4) 
r(136,4)=dyn(203,j) 
(137, })=dyn( 204, 4) 
r(138,4)=dyn(205, 4) 
r(308,j)=dyn(206,j) 
r(309,4)=dyn(207,j) 
r(139,4)=dyn( 208, j) 
r(140,4)=dyn(209, 4) 
(141, 4)=dyn(210,4) 
£(142,3)edyn( 211,34) 
£(310,4)=dyn( 212, 4) 
£(311,4)=dyn(213,4) 
r(143, j)=dyn(214,4) 
r(144,4)=dyn(215,4) 
£(145,4)=dyn(216,j) 
(146, 4)=dyn(217,3) 
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£(312,4)=dyn(218, 4) 
£(313, 3 )edyn(219, 4) 
r(147,4)=dyn(220,j) 
(148, 4)=dyn(221,4) 
r(149,4)«dyn(222,j) 
r(150,4)=dyn(223,3) 
£(314,4)=dyn(224, 4) 
r(315,4)=dyn(225,j) 
©(151, 4) =dyn( 226, 4) 
r(152,4)=dyn(227,j) 
r(153,4)edyn( 228, 4) 
£(154,4)=dyn(229, 4) 
r(316,j)=dyn(230,j) 
£(317,j)=dyn(231, 3) 
(155, j)=dyn(232,7) 
r(156,j)=dyn(233,3) 
£(157, 4) edyn(234, 3) 
r(158,4)=dyn(235,4) 
v(318,4)=dyn( 236, 7) 
r(319,j3)=dyn(237,3) 
6(159,4)=dyn(238,j3) 
r(160, j )=dyn( 239,35) 
r(161,j})=dyn(240, 3) 
(162, 4)=dyn(241,3) 
£(320,))=dyn(242,3) 
c(321,j)=dyn(243,j) 
©(163,4)=dyn(244,j) 
r(164, j)edyn(245, 4) 
r(165, 4) =dyn(246, 3) 
(166, 3)edyn( 247, 3) 
r(322,4)=dyn(248,3) 
(323, j)=dyn(249, 3) 
t(167,4)=dyn(250,3) 
(168, 4)=dyn(251,j) 
c(169,4)=dyn(252,3) 
r(170,4)=dyn(253,3) 
(324,43) =dyn(254,j) 
r(325,4)edyn( 255, 4) 
r(171, 3) =dyn(256, 4) 
(172, 4)=dyn( 257,74) 
(173, 4)=dyn(258,j) 
r(174,3)=dyn(259,j3) 
v(326,j)}=dyn( 260,34) 
(327, j)edyn({ 261, 4) 
(175, j)=dyn(262,j) 
r(176,5)=dyn(263,4) 
r(177, 4) =dyn(264,j) 
r(178,j})=dyn(265, 3) 
c(328,j4)=dyn(266, 4) 
r(329,4)=dyn( 267,34) 
r(179,4)=dyn( 268, j) 
c(180,j)=dyn( 269, }) 
c(181,4)=dyn(270,;) 
c(182,j)=dyn(271,j) 
£(330, 4}=dyn(272,3) 
(331, j)=dyn(273, 3) 
£(183,4)edyn(274,4) 
(184, 4)=dyn(275, 4) 
c(185,4)=dyn(276, 4) 
c(186,j)=dyn(277,3) 
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v(332,4)=dyn(278, 4) 
r(333,})=dyn(279, 4) 
c(187, 4)=dyn( 280,34) 
r(188, 4) =dyn( 281,34) 
(189, 4)=dyn( 282, 4) 
r(190,4)=dyn({283,j) 
£(334,4)=dyn( 284, 4) 
©(335,j)=dyn( 285, j) 
£(191,j)=dyn(286,j) 
£(192,4)=<dyn( 287, }) 
r(193,4)=dyn(288, 7) 
£(194,j)=dyn(289, 4) 
r(336,j)=dyn( 290, 4) 
£(337, 4)=dyn(291,3) 
£(195,4)=dyn(292,4) 
¢(196, 4) =dyn(293,j) 
£(197, j)=dyn(294,3) 
r(198,4)=dyn( 295, 4) 
£(338,4)=dyn(296,3) 
£(339,4)=dyn(297,3) 
r(199,4)=dyn( 298,34) 
£(200, j)=dyn(299,3) 
£(201, 3) =dyn( 300, 4) 
£(202,4)=dyn(301,j) 
(340, 4) =dyn( 302,34) 
©(341,4)=dyn( 303, 3) 
c(203,4)=dyn( 304,34) 
r(204,4)=dyn( 305, 3) 
£(205,4)=dyn( 306,35) 
©(206, j)=dyn( 307,34) 
0(342,4)=dyn( 308, j) 
c(343,4)=dyn(309, 3) 
r(207, j)=dyn( 310, 3) 
r(208,4)=dyn(311,3) 
r(209,4)=dyn(312,j) 
t(210,4)=dyn( 313, 3) 
(344, 5)=dyn(314,j) 
£(345,j)=dyn( 315,34) 
£(211,4)=dyn( 316,34) 
(212, 4)=dyn(317, 3) 
r(213,4)=dyn(318,j) 
r(214,3)=dyn(319,j4) 
(346, 4) =dyn(320,4) 
£(347,4)=dyn(321,3) 
r(215, j)=dyn( 322, 4) 
© (216, 4)=-dyn( 323,34) 
v(217,4)=dyn( 324, 3) 
£(218, 4)=dyn(325, 4) 
c(348,4)=dyn(326,j4) 
£(349,4)=dyn( 327,34) 
r(219,4)=dyn( 328, 4) 
r(220,4)=dyn(329,}4) 
t(221,4)=«dyn(330,j4) 
v(222,4)=dyn( 331,43) 
(350, 4) «dyn( 332,34) 
£(351,4)=dyn(333,4) 
(223, ))=dyn(334,3) 
t(224,4)=dyn( 335, 4) 
(225, 4)edyn( 336,35) 
c(226,4)=dyn( 337,35) 
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r(352,j)=dyn(338,4) 
£(353,4)=dyn(339,j) 
£(227,4)=dyn(340,j) 
(228,43) =dyn( 341,34) 
£(229,j)=dyn(342,j) 
£(230,4)=dyn(343,j) 
(354, j)=dyn(344,j) 
r(355, j)edyn(345,j) 
£(231,4)=dyn(346,3) 
r(232,j)=dyn(347,3) 
0(233, 4) =dyn( 348,35) 
(234, j)=dyn( 349, 4) 
c(356, 4) =dyn( 350, 4) 
£(357,4)=dyn( 351, 3) 
(235, j)=dyn(352,j3) 
£(236,j)=dyn(353,3) 
c(237, 4)=dyn(354,3) 
r(238,j)=dyn(355,j) 
r(358,4)=dyn(356,j) 
©(359,j)=dyn(357,j) 
©6(239,j)=dyn(358,j) 
(240, 4)=dyn( 359, j) 
r (360, 4) =dyn( 360, 3) 
r(361,j)=dyn( 361, 3) 
end do 


REARRANGE COLUMNS SO INTERNAL DISPLACEMENTS COME FIRST 


do i-1,361 
g(i,l)=r(i,1) 
g(i,2)=r(i,2) 


4)er(i,5) 


) 
44)er(i,14) 
45)er(i,15) 
1)=r(i,16) 
2)=c(i,17) 
3)er(i,18) 
714)er(i,19) 
1,246)ec(i,20) 
g(i,247)=r(i,21) 
g(i,15)er(i,22) 
g(i,16)=r(i,23) 
g(i,l7)=r(i,24) 
g(i,18)—c(i,25) 
g(1i,248)=r(i,26) 
g(i,249)=r(i,27) 
g(i,19)=r(i,28) 
g(i,20)=r(i,29) 
g(i,21)=c(i,30) 
g(i,22)=r(i,31) 
g(i,250)=r(i,32) 


woeunnnnonnennnne 
rte fade fade pute pate pate fate pate fede pute pete fete pete pete fete 
RERUN EP OOUNRaAU 
TY 
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g(i,251)=r(i,33) 
g(i,23)er(i,34) 
g(i,24)=r(i,35) 
g(i,25)er(i,36) 
g(i,26)er(i,37) 
g(i,252)er(i,38) 
g(i,253)er(i,39) 
g(i,27)=r(i,40) 
g(i,26)ec(i,41) 
gli, 29)er(i,42) 
g(i,30)=r(i,43) 
gli,254)ec(i,44) 
g(i,255)er(i,45) 
g(i,3l)=r(i, 46) 
g(i,32)=r(i,47) 
g(i,33)=c(i, 48) 
g(i,34)er(i,49) 
g(i,256)=r(i,50) 
g(i,257)er(i,51) 
g(i,35)=r(i,52) 
g(i,36)er(i,53) 
gli, 37)}=r(i,54) 
g(i,38)er(i,55) 
g(i,250)=r(i,56) 
g(i,259)=r(i,57) 
g(i,39)=c(i,58) 
g(i,40)er(i,59) 
g(i,41)=c(4,60) 
g(i,42)=r(i,61) 
g(i,260)=#r(i,62) 
g(i,261)=r(i,63) 
g(i,43)=rc(i,64) 
g(i,44)er(i,65) 
g(i,45)=c(i,66) 
g(i,46)=r(i,67) 
g(i,262)=r(i,68) 
g(i,263)=r(i,69) 
g(i,47)er(i,70) 
g(i,48)=r(i,71) 
g(i,49)er(i,72) 
g(i,50)=r(i,73) 
g(i,264)er(i,74) 
g(i,26S)er(i,75) 
g(i,51)=r(i,76) 
gli,52)=r(i,77) 
g(i,53)=r(i,78) 
g(i,54)=rc(i,79) 
g(i,266)=r(i,80) 
g(i,267)=r(i, 81) 
g(i,55)er(i,82) 
g(i,56)=r(i,83) 
g(i,57)=r(i,84) 
g(i,58)=c(i,85) 
g(i,268)=r(i,86) 
g(i,269)=r(i,87) 
g(i,59)=c(i,88) 
g(i,60)erc(i,89) 
g(i,61)#c(1,90) 
g(i,62)=c(i,91) 
g(i,270)er(i,92) 
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g(i,271)er(i,93) 
g(i,63)=r(4,94) 
g(i,64)=er(i,95) 
g(i,65)=r(i,96) 
g(i,66)er(i,97) 
g(i,272)=c(i,98) 
g(i,273)er(i,99) 
g(i,67)=£(i,100) 
g(i,68)=c(i,101) 
g(i,69)=c(i,102) 
g(i,70)=r(i,103) 
g(i,274)=r(i,104) 
g(i,275)=r(i,105) 
g(i,71)=£(i,106) 
g(i,72)=c(i,107) 
g(i,73)=c(i,108) 
g(i,74)=c(i,109) 
g(i,276)r(i,110) 
g(i,277)=c(i,111) 
g(i,75)ac(i,112) 
g(i,76)=c(i,113) 
g(i,77)=c(i,114) 
g(i,78)er(i,115) 
g(i,278)=r(i,116) 
g(i,279)=#r(i,117) 
g(i,79)ac(1,118) 
g(i,80)er(i,119) 
g(i,81)er(i,120) 
g(i,82)er(i,121) 
gli,280)er(i,122) 
g(i,281)er(i,123) 
g(i,83)=r(i,124) 
g(i,84)er(i,125) 
g(i,85)er(i,126) 
g(i,86)=r(i,127) 
g(i,282)=r(i,128) 
g(i,283)er(i,129) 
g(i,87)=r(i,130) 
g(i,88)=r(i,131) 
g(i,89)=r(i,132) 
g(i,90)=r(i,133) 
g(i,284)=r(i,134) 
g(i,285)=r(i,135) 
g(i,91)—r(i,136) 
g(i,92)=r(i,137) 
g(i,93)=r(i,138) 
g(i,94)=r(i,139) 
g(i,286)=r(i,140) 
g(i,287)=r(i,141) 
gli,95)er(i,142) 
gli,96)er(i,143) 
gli,97)er(i,144) 
gli, 98)ac(i,145) 
g(i,268)=r(i,146) 
g(i,289)=r(i1,147) 
g(i,99)=ar(i,148) 
g(i,100)=r(i,149) 
g(i,101)=c(i,150) 
g(i,102)er(i, 151) 
g(i,290)=r(i,152) 





41 











g(i,291)=r(i,153) 
g(i,103)=r(i,154) 
g(i,104)=r(i,155) 
g(i,105)er(i,156) 
g(i,106)=c(i,157) 
g(i,292)=r(i,158) 
gli, 293)=r(i,159) 
g(i,107)=c(i,160) 
g(i,208)=r(i,161) 
g(i,109)er(i,162) 
g(i,110)=r(i,163) 
g(i,294)=ac(i,164) 
g(i,295)=r(i,165) 
g(i,111)=c(i,166) 
g(i,122)"rc(i,167) 
g(i,113)=r(i,168) 
g(i,114)=r(i,169) 
g(i,296)=r(i,170) 
g(i,297)=#r(i,171) 
g(i,115)=r(i,172) 
g(i,116)=r(i,173) 
g(i,117)=r(i,174) 
g(i,118)=r(i,175) 
g(i,298)er(i,176) 
g(i,299)=c(i,177) 
g(i,119)=r(i,178) 
g(i,120)=r(i,179) 
g(i,121)=r(i,180) 
g(i,122)=c(i,181) 
g(i,300)er(i,182) 
g(i,301)er(i,183) 
g(i,123)=r(i,184) 
g(i,124)—r(i,185) 
g(i,125)=r(i,186) 
g(i,126)=r(i,187) 
g(i,302)=r(i,188) 
g(i,303)=c(i,189) 
g(i,127)=r(i,190) 
g(i,126)=r(i,191) 
g(i,129)=r(i,192) 
g(i,130)er(i,193) 
gli, 304)=r(i,194) 
g(i,305)=r(i,195) 
g(i,131)=c(i,196) 
g(i,132)er(i,197) 
g(i,133)=r(i,198) 
g(i,134)=r(i,199) 
g(i,306)=r(i,200) 
g(i,307)=r(i,201) 
g(i,135)er(i,202) 
g(i,136)=r(i,203) 
g(i,137)=r(i,204) 
g(i,138)er(i,205) 
g(i,308)=r(i,206) 
g(i,309)=r(i,207) 
g(i,139)=r(i,208) 
g(i,140)er(i,209) 
g(i,14l)=r(i,210) 
g(i,142)ec(i,211) 
g(i,310)er(i,212) 





g(i,31l)er(i,213) 
g(i,143)er(i,214) 
g(i,144)ar(i,215) 
g(i,145)er(i,216) 
g(i,146)=r(i,217) 
g(i,312)—r(i,218) 
g(i,313)—r(i,219) 
g(i,247)=r(i,220) 
g(i,148)=r(i,221) 
g(i,149)er(i,222) 
gii,150)=r(i,223) 
g(i,314)—c(i,224) 
g(i,315)=c(i,225) 
g(i,151)c(i,226) 
g(i,152)=c(i,227) 
g(i,153)=r(i,228) 
g(i,154)er(i,229) 
g(i, 316)=er(i,230) 
g(i,317)=r(i,231) 
g(i,155)=r(i, 232) 
g(i,156)=c(i,233) 
gli, 157 )=r(i,234) 
g(i,158)c(i,235) 
g(i,318)=r(i,236) 
g(i,319)=c(i,237) 
g(i,159)=c(i,238) 
g(i,160)=c(i,239) 
g(i,161)=c(i,240) 
g(i,162)er(i,241) 
g(i,320)=rc(i, 242) 
g(i,321)=r(i,243) 
g(i,163)ec(i,244) 
g(i,164)=r(i,245) 
g(i,165)=r(i,246) 
g(i,166)=r(i,247) 
g(i,322)—r(i,248) 
g(i,323)=r(i,249) 
g(i,167)—r(i,250) 
g(i,168)=r(i,251) 
g(i,169)=r(i,252) 
g(i,170)=r(i,253) 
g(i,324)er(i,254) 
g(i,325)-r(i,255) 
g(i,171)=r(i,256) 
g(i,172)er(i,257) 
g(i,173)=r(i,258) 
g(i,174)=r(i,259) 
g(i,326)=r(i, 260) 
g(i,327)—r(i,261) 
g(i,175)er(i,262) 
9(i,176)=r(i,263) 
g(i,177)=c(i, 264) 
g(i,178)=c(i,265) 
g(i,328)=r(i,266) 
g(i,329)=r(i, 267) 
g(i,179)=r(i,268) 
g(i,180)—r(i,269) 
g(i,181)=r(i,270) 
g(i,182)er(i,271) 
g(i,330)ec(i,272) 
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g(i,331)=r(i,273) 
g(i,183)=c(i,274) 
g(i,164)=c(i,275) 
g(i,185)=r(i,276) 
g(i,186)=r(i,277) S 
g(i,332)er(i,278) 
g(i,333)=c(i,279) 
g(i,187)=r(i, 280) 
g(i,168)=r(i,281) . 
g(i,189)=r(i,282) 
g(i,190)=c(i, 283) 
g(i,334)=er(i, 284) 
g(i,335)=c(i,285) 
g(i,191)er(i,286) 
g(i,192)=r(i, 287) 
g(i,193)=r(i, 288) 
g(i,194)=r(1,289) 
g(i,336)=r(i,290) 
g(i,337)er(i,291) 
g(i,195)-r(i,292) 
g(i,196)=r(i,293) 
g(i,197)=r(i,294) 
9(i,198)=r(i, 295) 
g(i,338)=r(i,296) 
g(i,339)=c(i,297) 
g(i,199)=r(i,298) 
g(i,200)=r(i,299) 
g(i,201)=r(i, 300) 
g(i,202)=r(i, 301) 
g(i,340)=c(i, 302) 
g(i,341)=#r(i, 303) 
g(i,203)=r(i, 304) 
g(i,204)=c(i,305) 
g(i,205)=r(i, 306) é 
g(i,206)=r(i, 307) 
g(i,342)=c(i, 308) 
g(i,343)=r(i, 309) 
g(i, 207 )ac(i, 310) 
g(i,208)=r(i, 311) - 
g(i,209)r(i, 312) 
g(i,210)er(i, 313) 
g(i,344)er(i,314) 
g(i,345)=r(i,315) 
g(i,211)=r(i,316) 
g(i,212)=c(i, 317) 
g(i,213)er(i,316) 
g(i,214)=r(i,319) 
g(i,346)=c(i, 320) 
g(i,347)ac(i,321) 
g(i,215)=c(i, 322) 
g(i,216)=r(i,323) 
g(i,217)=r(i, 324) 
g(i,218)=r(i, 325) 
g(i,348)=c(i, 326) 
gli, 349)=c(i,327) 
g(i,219)er(i, 328) 
g(i,220)=r(i, 329) 
g(i,221)—r(i,330) 
g(i,222)=r(i,331) 
gi, 350)ec(i, 332) 
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g(i,351)=c(i,333) 
g(i,223)=—r(i,334) 
g(i,224)=r(i,335) 
g(i,225)=c(i, 336) 
g(i,226)=r(i,337) 
g(i,352)=r(i, 338) 
g(i,353)=r(i,339) 
g(i,227)=c(i, 340) 
g(i,228)er(i,341) 
g(i,229)=r(i,342) 
g(i,230)=c(i,343) 
g(i,354)=r(i, 344) 
g(i,355)er(i, 345) 
g(i,231)=r(i,346) 
g(i,232)—r(i,347) 
g(i,233)er(i,348) 
g(i,234)=r(i, 349) 
g(i,356)=r(i, 350) 
g(i,357)=r(i,351) 
g(i,235)er(i, 352) 
9(i,236)er(i, 353) 
g(i,237)=r(i,354) 
g(i,238)=r(i,355) 
g(i,358)=c(i,356) 
g(i,359)=r(i, 357) 
g(i,239)=r(i,358) 
g(i,240)=r(i, 359) 
g(i,360)=r(i, 360) 
g(i,361)=r(i,361) 
end do 


ASSEMBLE TRANSFORMATION MATRIX, t( 361,361) 
WHICH CONVERTS SURFACE DOFs TO NORMAL AND TANGENTIAL 


Upper left corner 
do i=1,240 
do j=1,240 
t(i,j)=0.0 
end do 
t(i,i)=1.0 
end do 


Upper right corner 
do i=1,240 
do j=241,361 
t(i,j)=0.0 
end do 
end do 


Lower left corner 
do i=241,361 
do j=1,240 
t(i,j)=0.0 
end do 
end do 


Lower cight corner 
do i=241,361 
do j=-241,361 
t(i,j)=0.0 
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end do 
end do 
t(241,241}=-1.0 
t(360,360)=1.0 
t(361,361)=1.0 


READ IN SURFACE NODE COORDINATES AND COMPUTE SIN THETA AND COS THETA 
FOR ALL BUT FIRST AND LAST SURFACE NODES 
NOTE THETA IS THE POLAR ANGLE 
do i=1,59 
read(5,°(2£10.8)'’) xx,yy 
trmesqrt(xx**2+yy**2) 
ssin(i)=xx/rr 
ccos(i)=yy/rr 
end do 


k=241 

do 1=242,358,2 
t(i,i)#ssin(i-k) 
t(i,i+1)=<-ccos(i-k) 
t(i+1,i)=<ccos(i-k) 
t(i+l,i+l)=ssin(i-k) 
kek+1 

end do 


ASSEMBLE tt(361,361) = TRANSPOSE OF t(361, 361) 
do i=1,361 
do j=1,361 
tt(j,i)=t(i,j) 
end do 
end do 


MULTIPLY g BY t 
do i=1,361 
do j~1,361 
sum=0.0 
do k=1,361 
sum=sumtg(i,k)*t(k,j) 
end do 
pint(i,j)=sum 
end do 
end do 


do i=1,361 
do j=1,361 
sum=0.0 
do k=1,361 
sum=sum+tt(i,k)*pint(k,j) 
end do 
prod(i,j)=sum 
end do 
end do 


TEST SYMMETRY OF TRANSFORMED FULL DYNAMICAL MATRIX 
nerr=0 
do i-1,361 
do j-i,361 
if (prod(j,i).ne.prod(i,j)) then 
write (*,*)*Symm error in transf full dyn mtrx at (i,j)#'’,i,j 
nerrenerr+l 
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end if 
end do 
end do 
write (*,*) ‘Number of symmetry errors in transf full dyn mtrx ='’,nerr 
4£ (merr.gt.0) then 
STOP 
end if 


REARRANGE ROWS SO ZERO (TANGENTIAL) SURFACE FORCES COME FIRST 
do i=1,240 

do j=1,361 

s(i,j)=prod(i,j) 

end do 
end do 
do j=1,361 
s(300,j)=prod(241,j) 
8(301,j)=prod(242,j) 
6(241,j})=prod(243, 3) 
s(302,j)=prod(244,j) 
8(242,j)=prod(245,j) 
s(303,})=prod(246,j) 
s(243,j})=prod(247, 3) 
s (304, j)=prod(248, j) 
s(244,j)=prod(249,j) 
s(305,4)=prod(250,j) 
s(245,j})=prod(251,j) 
s(306,j)=prod(252,j) 
s(246,j)=prod(253,j) 
s(307,3)=prod(254,3) 
8(247,j)=prod(255,j) 
s(308,j)=prod(256, 3) 
s(248,j})=prod(257,j) 
s(309,j)=prod(258,j) 
8(249,j)=prod(259,j) 
8(310,j)=prod(260,j) 
s(250,4)=prod( 261, 3} 
8(311,j)=prod(262,j) 
s(251,j)=prod(263,j) 
s(312,j)=prod(264,j) 
s(252,j)=prod(265,}) 
s(313,4)=prod(266,j) 
s(253,j)=prod(267,j) 
s(314,4)=<prod(268,j) 
§(254,j)=prod(269,j) 
s(315,j)=prod(270,j) 
s(255,j)=prod(271,j) 
s(316,})=prod(272,3) 
8(256,4)=prod(273,j) 
s(317,4)=prod(274,}) 
s(257,})=pcrod( 275, 4) 
s(3186,j)=prod(276,j) 
8(258,j)=prod(277,j) 
s(319,})=prod(278, 3) 
8(259,4)=prod(279,j) 
8(320,4)=prod( 260, 4) 
§(260,4)=prod(261,j) 
8(321,j4)=prod( 282, 4) 
8(261,4)=prod( 283, 3) 
s(322,j)=prod(284,j) 
s(262,4)=prod(285,j) 
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s(323,4))=prod(286,j 
8(263,})=prod(287,j 
8(324,})=prod(268,j 
8(264, })=prod(289,4 
(325, 4)=prod(290,j 
s(265,j)=prod(291,j 
5(326, })=prod(292,j 
s(266,4)=prod(293,j 
s(327,4j)=prod(294,4 
s(267,4)=prod(295,j 
8(328, 3) =prod(296,3 
6(268,j)=prod(297,j 
s(329,})=prod(298,j 
s(269,j)=prod(299,j 
s(330, j)=prod( 300, 3 
s(270,4)=prod(301,j 
s(331,4)=prod(302,j 
8(271, 4) =prod(303,j 
s(332,4)=prod(304,j 
s(272,}4)=prod(305,j 
s(333,4)=prod(306,j 
s(273, j)=prod( 307, 
s(334,4)=prod(308,j 
(274, j)=prod{ 309, j 
s(335,4)=prod(310,j 
s(275, })=prod(311,j 
s (336, })=prod( 312, 
5 (276, j)=prod(313,j 
(337, j})=prod(314,j 
s(277,4)=prod(315,j 
s(338,j)=prod(316,j 
s(278,j})=prod(317,j 
8(339, j)=prod(318,j 
8(279,j)=prod( 319, 4 
6(340, j)=prod(320,j 
6(280, 3) =prod( 321, 3 
6(341,})=prod(322,j 
s(281. =prod(323,j 
s(342,3, =prod(324,j 
s(282, j)=prod(325,j 
8(343,j)=prod(326,j 
6(283,})=prod(327,j 
$(344, j)=prod( 328,j 
5(284, j)=prod(329,j 
6(345,4)=prod(330,) 
s(285, j)=pcod(331,j 
6(346, })=prod(332,3 
8(286, 4) =prod(333,j 
8(347,j)=prod(334,j 
6(287,4)=prod(335,j 
6(348,})=prod( 336,j 
s(288, })=prod( 337, 
8( 349, j)=prod(338,3 
6(289,})=prod(339,j 
8(350,j)=prod( 340, 3 
8(290,j)=prod( 341,j 
6(351,4)=prod(342,3 
8(291,4)=prod(343,j 
6(352,j)=prod( 344,j 
8(292,j)=prod( 345, 


hate tedetdete tee eth ee eee 


mee eee ee eae ee et et ee ee ee ee ee 
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8(353,j)=prod( 346, j 
8(293,j)=prod( 347, j 
8(354,j)=prod(348,j 
8(294,))=prod(349, 3 
8(355,4)=prod(350,j 
8(295,j})=prod(351,j 
s(356,4)=prod(352,j 
8(296,j)=prod(353,j 
8(357,4)=prod(354,j 
8(297,4)=prod(355, j 
8(358,j)=prod(356,j 
s(298,4)=prod(357,j 
8(359,4)=prod(358,j 
8(299,4)=prod(359,j 
8(360,4)=prod(360,j 
s(361,4)=prod(361,j 
end do 


REARRANGE COLUMNS SO TANGENTIAL SURFACE DISPLACEMENTS COME FIRST 


do i=1,361 
do j=1,240 


end d 
end do 
do i=1,361 
mat(i,300)=s(i,241) 
mat(i,301)<#s(i,242) 
mat(i,241)=#s(i,243) 
mat(i,302)«s(i,244) 
mat(i,242)=s(i,245) 
mat(i,303)=s(i,246) 
mat(i,243)=s(i,247) 
mat(i,304)=s(i,248) 
mat(i,244)=s(i,249) 
mat(i,305)=s(i,250) 
mat(i,245)=#s(i,251) 
mat(i,306)=s(i,252) 
mat(i,246)=s(i,253) 
mat(i,307)=s(i,254) 
mat(i,247)=s(i,255) 
mat(i,308)=s(i,256) 
mat(i,248)=#s(i,257) 
mat(i,309)=s(i,258) 
mat(i,249)=s(i,259) 
mat(i,310)=s(i,260) 
mat(i,250)=s(i,261) 
mat(i,311)#s(i,262) 
mat(i,251)=s(i,263) 
mat(i,312)<s(i,264) 
mat(i,252)=s(i,265) 
mat(i,313)=s(i,266) 
mat(i,253)=s(i,267) 
mat(i,314)=s(i,268) 
mat(i,254)es(i,269) 
mat(i,315)es(1,270) 
mat(i,255)=s(i,271) 
mat(i,316)=s(i,272) 
mat(i,256)=#s(i,273) 
mat(i,317)s(i,274) 
mat(i,257)=s(i,275) 


) 
) 
) 
) 
) 
) 
) 
) 
) 
) 
) 
) 
) 
) 
) 
) 


mat(i,j)<s(i,j) 
° 
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mat(i,316)=s(i,276) 
mat(i,258)=s(i,277) 
mat(i,319)-s(i,278) 
mat(i,259)=s(i,279) 
mat(i,320)=s(i,280) 
mat(i,260)-s(i,281) 
mat(i,321)=s(i,262) 
mat(i,261)=s(i,283) 
mat(i,322)=<s(i,284) 
mat(i,262)=s(i,285) 
mat(i,323)=s(i,286) 
mat(i,263)=s(i,287) 
mat(i,324)<s(i,268) 
mat(i,264)=s(i,289) 
mat(i,325)=s(i,290) 
mat(i,265)=s(i,291) 
mat({i,326)=s(i,292) 
mat(i,266)=s(i,293) 
mat(i,327)=s(i,294) 
mat(i,267)=s(i,295) 
mat(i,328)=s(i,296) 
mat(i,268)=s(i,297) 
mat(i,329)=s(i,298) 
mat(i,269)=s(i,299) 
mat(i,330)=s(i,300) 
mat(i,270)=s(i,301) 
mat(i,331)=s(i,302) 
mat(i,271)=<s(i,303) 
mat(i,332)=s(i,304) 
mat(i,272)=#s(i,305) 
mat(i,333)=s(i,306) 
mat(i,273)=s(i,307) 
mat(i,334)=s(i, 308) 
mat(i,274)=s(i,309) 
mat(i,335)=s(i,310) 
mat(i,275)es(i,311) 
mat(i,336)=s(i,312) 
mat(i,276)=s(i,313) 
mat(i,337)=s(i,314) 
mat(i,277)=s(i,315) 
mat(i,338)=s(i,316) 
mat(i,278)=<s(i,317) 
mat(i,339)=s(i,318) 
mat(i,279)=s(i,319) 
mat(i,340)=s(i,320) 
mat(i,280)=s(i,321) 
mat(i,342)=2s(i,322) 
mat(i,281)=s(i,323) 
mat(i,342)=s(i,324) 
mat(i,282)=<s(i,325) 
mat(i,343)=s(i,326) 
mat(i,283)=8s(i,327) 
mat(i,344)«s(1i,328) 
mat(i,284)-s(i,329)} 
mat(i,345)=s(i,330) 
mat(i,285)=s(1i,331) 
mat(i,346)=s(i,332) 
mat(i,286)=s(4i,333) 
mat(i,347)es(i,334) 
mat(i,287)=s(i,335) 
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mat(i,346)=s(i,336) 
mat(i,288)es(i,337) 
mat(i,349)=s(i,338) 
mat(i,289)=s(i,339) 
mat(i,350)=s(1,340) 
mat(i,290)=s(i,341) 
mat(i,351)«s(i,342) 
mat(i,291)=s(i,343) 
mat(i,352)=s(i,344) 
mat(i,292)=s(i,345) 
mat(i,353)=s(i,346) 
mat(i,293)=s(i,347) 
mat(i,354)<s(i,348) 
mat(i,294)=s(i,349) 
mat(i,355)=s(i,350) 
mat(i,295)=#s(i,351) 
mat(i,356)=s(i,352) 
mat(i,296)=s(i,353) 
mat(i,357)=s(i,354) 
mat(i,297)=s(i,355) 
mat(i,358)=s(i,356) 
mat(i,298)=s(i,357) 
mat(i,359)=s(i,358) 
mat(i,299)=s(i,359) 
mat(i, 360)=s(i, 360) 
mat(i,361)<s(i,361) 
end do 


PARTITION mat MATRIX INTO a,b,c AND d MATRICES 
do i#1,299 
do j=1,299 
a(i,j)=mat(i,j) 
end do 
end do 


do i=1,299 
do j=300, 361 
b(i,j-299)=mat(i,j) 
end do 
end do 


do i=300,361 
do j=1,299 
c(i-299,4)=mat(i,j) 
end do 
end do 


do i=300,361 
do j=300,361 
d(i-299, j~299)=mat(i,j) 
end do 
end do 


COMPUTE INVERSE OF MATRIX a 
n=299 
np=299 
do i=l,n 
do jel1,n 
y(i,j)=0.0 
end do 
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y(i,i)=1.0 
end do 
call ludcmp(a,n,np, indx,d) 
do j=i,n 
call lubksb(a,n,np,indx,y(1,4)) 
end do 


MULTIPLY c BY y 
do i=1,62 
do j#1,299 
sum=0.0 
do k=1,299 
sum=sumtc(i,k)*y(k,j) 
end do 
pli.j)=sum 
end do 
end do 


MULTIPLY p BY b 
do i=1,62 
do j=1,62 
sum=0.0 
do ke1,299 
sum=sum+p(i,k)*b(k, 4) 
end do 
pp(i,j)=sum 
end do 
end do 


COMPUTE d-pp 
do i=1,62 
do j=1,62 
nn(i,j)=d(i,j)-ppli,j) 
end do 
end do 


TEST SYMMETRY OF REDUCED DYNAMICAL MATRIX 
nerr=0 
do i=1,62 
do j=i,62 
fracerr=Abs(nn(j,i)-nn(i,j))/Sqrt(Abs(nn(j,i)*nn(i,j))) 
if (fracerr.ge.0.001) then 
write (*,*)’Symmetry error in reduced dyn. mtrx at (i,j)=',i,j 
necc=nerrctl 
end if 
end do 
end do 
write (*,*) "Number of symmetry errors in reduced dyn. matrix ='’,nerr 
if (nerr.gt.0) then 
STOP 
end if 


THIS COMPLETES THE REDUCTION TO NODAL SURF VEL DEG OF FREEDOM 
STORE REDUCED DYNAMICAL MATRIX IN UNFORMATTED BINARY FORM 
do i-1,62 
weite(4)(nn(i,j),j=<1,62) 
end do 


end 
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END OF MAIN PROGRAM 


subroutine ludcmp(a,n,np, indx,d) 
parameter (nmax=299, tiny=1 .0e-32) 
dimension a(299,299), indx(299), vv(299) 
d=1.0 
do i-l1,n 
aamax=0.0 
do j=i,n 
if(abs(a(i,j)).gt.aamax) aamaxeabs(a(i,j)) 
end do 
if(aamax.eq.0.0) pause ‘singular matrix.’ 
vv(i)=1.0/aamax 
end do 


do j=l,n 

do i=-1,j-1 
sum=a(i,j) 
do ke#1,i-1 

sum=sum-a(i,k)*a(k,4) 

end do 
a(i,j)=sum 

end do 

aamax=0.0 


do i=j,n 
sum~a(i,j) 
do k#1,j-1 
sum=sum-a(i,k)*a(k,j) 
end do 
a(i,j)<sum 
dum=vv(i)*abs (sum) 
if(dum.ge.aamax) then 
imax=i 
aamax=dum 
end if 
end do 


if(j.ne.imax)then 
do k=1,n 
dum=a(imax,k) 
a(imax,k)=a(4,k) 
a(j,k)=dum 
end do 
d=-d 
vv(imax)=vv(4) 
end if 
indx(j)=imax 
if(a(j,j)-eq.0.) a(j,j)=tiny 
if(j.ne.n) then 


dum=1.0/a(j,j) 
do i#j+l,n 
a(i,j)ea(i,j)*dum 
end do 
end if 

end do 

return 

end 
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subroutine lubksb(a,n,np,indx,hi) 
cimenetee a(299,299), indx( 299) ,hi(299) 
ii= 
do ii,n 

m=indx(i) 

sum-hi(m) 

hi(m)=hi(i) 

if(ii.ne.0) then 

do jeii,i-1 

sum=sum-a(i,j)*hi(j) 

end do 

else if(sum.ne.0.) then 

iisi 

end if 

hi(i)=sum 
end do 


do ien,1,-1 
sum=hi(i) 
if(i.lt.n) then 
do j=i+i,n 
sum=sum-a(i,j)*hi(j) 
end do 
end if 
hi(i)=sum/a(i,i) 
end do 
ceturn 
end 
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APPENDIX C: 


PROGRAM SPHMTRX: rads in reduced nodal dynamical matrix and produces 
dynamical matrix in terms of first seven Legendre polynomials 


real pi,sum 
real xx,yy,rr,ccos(61) 
real Dnod(62,62),P(61,7),Cinv(7,7),Pt(7,61),Dsph(8,8) 


OLD FILES 
open(unit=8,file=’transform.sph’ ,status='old’ ) 
open(unit=9,file="reddynmtrx.dat’,status#’old’,form=’unformatted’ ) 


NEW FILES 

open(unit=12,file='Dsph.dat’,status='new’ ,form="unformatted’ ) 
open(unit=13,file="Dsphl.lis’,status='new’ ) 
open(unit=14,file='Dsph2.lis’,status=’new’? ) 


pi=acos(-1.0) 


READ IN COORDINATES OF ALL EXCEPT FIRST AND LAST SURFACE NODES 
AND COMPUTE COS OF EACH ANGLE. 
Note first node is at theta-pi and last node is at theta=-0 where 
theta ia the polar angle 
do i#1,59 
read(8,‘(2£10.8)') xx,yy 
trmsqrt(xx**2+yy**2) 
ccos(itl)eyy/rcr 
end do 
ccos(1)=#-1.0 
ccos(61)=1.0 


UNFORMATTED READ IN REDUCED DYNAMICAL MATRIX Dnod(62,62) 
do i=1,62 

read(9)(Dnod(i,j),j=1,62) 
end do 


COMPUTE NODAL DISPLACEMENT TRANSFORMATION MATRIX, P(61,7), USING THE 
FIRST 7 LEGENDRE POLYNOMIALS EVALUATED AT COS(NODE3) THRU COS(NODE163) 
Compute zeroth through sixth order Legendre polynomials 
do i=1,61 

c2=ccos(i)**2 

P(i,1)=1.0 

P(i,2)=ccos(i) 


P(i,3)=(3.0*c2-1.0)/2.0 
P(i,4)=((5.0*%c2-3.0)*ccos(i))}/2.0 
P(i,5)=((35.0%c2~30.0)*c2+3.0)/8.0 
P(i,6)=(((63.0%c2-70.0)*c2+15.0)*ccos(i))/8.0 
Oa a aa el al ea ae oa 
end do 


COMPUTE TRANSPOSE OF NODAL DISPLACEMENT TRANSFORMATION MATRIX Pt(7,61) 
do i#1,7 
do j=1,61 
Pt(i,j)=P(4j,1) 
end do 
end do 


COMPUTE INVERSE OF LEGENDRE POLYNOMIAL NORMALIZATION CONSTANTS Cinv(7,7) 
do i-1,7 
do j#1,7 
Cinv(i,j)=#0.0 
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end do 
Cinv(i,i)=2.0/(2*(i-1)41) 
end do 


tints UPPER LEFT 7x7 BLOCK OF Dsph(8,8)=Cinv(7,7)*Pt(7,62)*Dnod(61,61) ~ 
o i@l,7 I 
do j=1,7 ; 
sum=0.0 
do k=1,7 
do m=1,61 7 
do n«1,61 
sum=sum+Cinv(i,k)*Pt(k,m)*Dnod(m,n)*P(n,j) 
end do 
end do 
end do 
Dsph(i,j)=sum 
end do 
end do 





Aarts UPPER RIGHT 7x1 BLOCK OF Dsph(8,8)=Cinv(7,7)*Pt(7,61)*knod(61,1) 
do i#1,7 
sum=0.0 
do j=1,7 
do k«1,61 
sum=sum+Cinv(i,j)*Pt( j,k) *Dnod(k,62) 
end do 
end do 
Dsph(i,8)=sum 
end do 


COMPUTE LOWER LEFT 1x? BLOCK OF Dsph(8,8)=Knodt(1,61)*P(61,7) 
do j=1,7 

sum=0.0 

do k#1,61 

sum=sum+Dnod(62,k)*P(k,j) 

end do 

Dsph(8,j)=sum 
end do 


BLOCKED CAPACITANCE IS UNCHANGED 
Dsph(8,8)=Dnod(62,62) 


STORE Dsph(8,8) 


do iwi,8 
write(12)(Dsph(i,j),j=1,8) 
write(13,10)(Dsph(i,j),je1,4) 
write(14,10)(Dsph(i,j),j=5,8) 

end do 

format(3x,4e17.7,//) 


end 
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